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ABSTRACT 

Context. Radiative feedback plays a crucial role in the formation of massive stars. The implementation of a fast and accurate 
description of the proceeding thermodynamics in pre-stellar cores and evolving accretion disks is therefore a main effort in current 
hydrodynamics simulations. 

Aims. We introduce our newly implemented three-dimensional frequency dependent radiation transport algorithm for hydrodynamics 
simulations of spatial configurations with a dominant central source. 

Methods. The module combines the advantage of the speed of an approximate Flux Limited Diffusion (FLD) solver in the one- 
temperature approach, which is valid in the static diffusion limit, with the high accuracy of a frequency dependent first order ray- 
tracing routine. The ray-tracing routine especially compensates the introduced inaccuracies by standard approximate FLD solvers in 
transition regions from optically thin to thick and yields the correct optical depths for the frequency dependent stellar irradiation. 
Both components of our module make use of realistic tabulated dust opacities. The module is parallelized for distributed memory 
machines based on the message passing interface standard. We implemented the module in the three-dimensional high-order magneto- 
hydrodynamics code Pluto. 

Results. We prove the viability of the scheme in a standard radiation benchmark test compared to a full frequency dependent Monte- 
Carlo based radiative transfer code. The setup includes a central star, a circumstellar flared disk, as well as an envelope. The test is 
performed for different optical depths. Considering the frequency dependence of the stellar irradiation, the temperature distributions 
can be described precisely in the optically thin, thick, and irradiated transition regions. Resulting radiative forces onto dust grains are 
reproduced with high accuracy. The achievable parallel speedup of the method imposes no restriction on further radiative (magneto-) 
hydrodynamics simulations. 

Conclusions. The proposed approximate radiation transport method enables frequency dependent radiation hydrodynamics studies 
of the evolution of pre-stellar cores and circumstellar accretion disks around an evolving massive star in a highly efficient and accurate 
manner. 
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1. Introduction 

Massive stars are of importance for a wide range of astrophys- 
ical problems. In spite of the fact that they are both rare and 
short-living, they represent the major source of radiation energy 
in their stellar clusters. Therefore, they act as valuable tracers 
of star formation rates in distant galaxies. Additionally, the ra- 
diation emitted during a massive star's lifetime influences the 
surroundings through various interactions such as heating and 
ionization of gas, evaporation of dust, and radiative forces lead- 
ing to powerful stellar winds and outflows. Finally, when a mas- 
sive star dies, it enriches its neighborhood with heavy elements. 
In this sense, massive stars are the main drivers of the morpho- 
logical, dynamical, and chemical evolution of their complex en- 
vironments (Zinnecker & Yorke 2007; McKee & Ostriker 2007; 
Henning 2008). 

However, knowledge about their formation is rather 
poor compared to the case of low-mass star formation. 
Observationally, this is mostly due to their larger average dis- 
tances and the fact that they are deeply embedded in dense, 



opaque cores especially during their early evolutionary phases. 
Also, their short lifetime, rareness, and complex environment 
pose difficulties for detailed observations. Nevertheless, cur- 
rent observational results, e.g. from the Spitzer Space Telescope 
Survey and millimeter interferometry, support the assumption 
that the basic concepts of star formation, including the collapse 
of an unstable gas and dust core, the formation of both bipolar 
outflows and jets, as well as accretion disks, are possibly valid 
throughout the whole range of stellar masses (e.g. Bally 2008; 
Bik et al. 2008; Kumar & Grave 2008; Steinacker et al. 2008; 
VehofF et al. 2008). A good review of observational results can 
be found for instance in Beuther et al. (2007). Future genera- 
tions of space telescopes and interferometric systems such as the 
Herschel Space Observatory, the James Webb Space Telescope 
(JWST), and the Atacama Large Millimeter Array (ALMA) will 
provide a deeper insight into the mechanisms of in-falling en- 
velopes, bipolar outflows, and massive accretion disks. This will 
definitely put tighter constraints on current theoretical models. 
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From the theoretical point of view, assuming that the for- 
mation of high-mass stars is basically a scaled-up version of 
low-mass star formation, implies its own challenges: Their rapid 
evolution, especially the shorter Kelvin-Helmholtz contraction 
timescale, leads to the interaction of the accreting flow of gas 
and dust with the emitted radiation from the newborn star (Shu 
et al. 1987). In the extreme case of a purely spherically sym- 
metric collapse of a pre-stellar core (without angular momen- 
tum), the in-fall is potentially stopped by the growing radiation 
pressure entirely. In the static limit, the equilibrium between 
radiative and gravitational forces at the dust condensation ra- 
dius is described by the Eddington limit L*/M* = (4n G c)/k,, 
where L„, M», and k* denote the stellar luminosity, mass, and 
the dust opacity regarding the stellar irradiation spectrum, re- 
spectively. The underlying one-dimensional collapse problem is 
extensively studied in analytical studies (Penston 1969b; Shu 
1977; Hunter 1977) as well as in numerical simulations (Penston 
1969a; Hunter 1977; Whitworfh & Summers 1985; Foster & 
Chevalier 1993), either in hydrodynamics or as evolution of self- 
similar solutions. Also stellar radiative feedback has been con- 
sidered in one-dimensional studies (Larson & Starrfield 1971; 
Kahn 1974; Yorke & Kriigel 1977; Wolfire & Cassinelli 1987). 
However star formation is rarely a perfectly spherically symmet- 
ric process. 

As stated by Nakano (1989), massive star formation is at 
least a two-dimensional non-spherically symmetric problem, in- 
cluding initial angular momentum to force the formation of a 
circumstellar disk and polar cavities. Both will potentially help 
to overcome the Eddington barrier. On small scales, Turner et al. 
(2007) show that the so-called photon bubble instability can sig- 
nificantly enhance accretion and avoid the radiation pressure 
problem. On large scales, the formation of bipolar cavities pos- 
sibly reduces the radiative force onto the accretion flow by the 
so-called flashlight effect (Yorke & Sonnhalter 2002; Krumholz 
et al. 2005b). Although Krumholz et al. (2005a, 2007a, 2009) 
claim that the Eddington barrier may be broken via a Rayleigh- 
Taylor instability, it is unclear if this instability is the most im- 
portant one and how this instability will change in case of real- 
istic frequency dependent radiative feedback instead of the used 
gray (i.e. frequency averaged) Flux Limited Diffusion (hereafter 
called FLD) approximation. 

The important role of frequency dependent radiation trans- 
port in the formation of massive stars was already shown in 
Yorke & Sonnhalter (2002), but due to the huge computational 
overhead of the frequency dependent radiation transfer it was 
neither possible to study a huge amount of different initial con- 
ditions (to scan the parameter space), nor to perform high- 
resolution simulations of the accretion process. 

Regarding these issues, we demonstrate here the usability of 
our newly developed fast three-dimensional frequency depen- 
dent approximate radiation transport module for numerical hy- 
drodynamics. 

The most accurate description of the physics proceeding 
during the collapse of a pre-stellar core would include a fre- 
quency dependent radiation transport step after each hydro- 
dynamic timestep using a modern Monte-Carlo based or ray- 
tracing radiative transfer method. The CPU time needed to solve 
one single hydrodynamic timestep is generally orders of mag- 
nitude lower than the time for a frequency dependent radiation 
transport step, especially in complex geometries. Thus, this ap- 
proach is not feasible with current computing technology for a 
large number of grid cells in more than one dimension due to 
the large amount of computational time needed for each radi- 
ation transport step. A more desired radiation hydrodynamics 



scheme should roughly spend the same CPU time on the radia- 
tion physics as on the hydrodynamics part. Sensible approxima- 
tions are thus necessary to speed up the radiative transfer in such 
hydrodynamics studies. 

In hydrostatic disk atmosphere computations, Murray et al. 
(1994) introduced a splitting of a two-dimensional radiation field 
into an irradiated and a diffuse part. Already Wolfire & Cassinelli 
(1986, 1987) used such a splitting in one dimension to study the 
radiation feedback of massive stars, which has been shown to be 
a valid approach in one-dimensional hydrodynamics simulations 
by Edgar & Clarke (2003). In the present paper, we expand the 
method to higher dimensions and show its validity in an axially 
symmetric setup, consisting of a central star, a flared circumstel- 
lar disk, and an envelope. In fact, we found that it is necessary 
to consider the frequency dependence of the stellar irradiation 
feedback to reconstruct a reasonable approximation to the spa- 
tial temperature distribution. 

A first approach for fast two-dimensional axially symmet- 
ric radiative transfer is, for example, the gray diffusion approx- 
imation studied by Tscharnuter & Boss (1993), which is only 
applicable in the optically thick limit. Ray-tracing based meth- 
ods (like in Efstathiou & Rowan-Robinson 1990) show high ac- 
curacy, but as already mentioned also require much CPU time, 
which yields low efficiency in combination with hydrodynamics 
solvers. Another common approach for the description of radia- 
tive processes in circumstellar disks is the gray FLD approxi- 
mation (Kley 1989; Bodenheimer et al. 1990; Klahr et al. 1999; 
Klahr & Kley 2006). It provides a fast method to determine the 
temperature evolution in the optically thick (diffusion) as well 
as in the optically thin (free-floating) limit, but shows stronger 
deviations in transition regions, demonstrated e.g. in Boley et al. 
(2007). In the case of gray FLD, which is still the default tech- 
nique in modern radiation hydrodynamics codes, this method of 
course suffers strongly from the lack of frequency dependence, 
when compared to the accuracy of modern Monte-Carlo based 
codes. 

In the following, we present our results in combining the 
advantage of the gray FLD approximation (speed) with the 
accuracy of frequency dependent ray-tracing. The approxima- 
tion described here results in a large reduction in computing 
time compared to Monte-Carlo based radiative transfer. This 
allows us to implement this particular code in the framework 
of three-dimensional (magneto-) hydrodynamics simulations of 
circumstellar disks and in-falling envelopes on a parallel decom- 
posed (spherical) grid. In particular, we implemented our ra- 
diation transport module into the three-dimensional high-order 
Godunov magneto-hydrodynamics code Pluto (Mignone et al. 
2007). This tool allows us to perform high-resolution studies of 
accretion flows around a forming massive star. 

The basic equations of the approximate radiation transport 
method and the numerical details are described in the follow- 
ing Sect. 2. The test of the accuracy of the proposed radia- 
tion transport scheme in the setup of a currently standard ra- 
diation benchmark test (Pascucci et al. 2004) is presented in 
Sect. 3. The parallel scaling behavior of our implementation is 
discussed in Sect. 4. In Sect. 5, we show the results of two stan- 
dard radiation hydrodynamics shock tests performed with the 
FLD part of our implemented radiation transport method and 
the open source magneto-hydrodynamics code Pluto (Mignone 
et al. 2007). Conclusions as well as further studies planned are 
described in Sect. 6. 
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2. Theory and numerics of the approximate 
radiation transport scheme 

In the following, we recapitulate the general ideas, the basic 
equations, and methods of the frequency dependent approxi- 
mate radiation transport scheme. This section should allow the 
reader to follow our motivation for the newly implemented three- 
dimensional radiative transfer module. Every implemented for- 
mula and numerical detail is given. 

The general idea of the method is to split the radiation trans- 
port into two components (Wolfire & Cassinelli 1986; Murray 
et al. 1994; Edgar & Clarke 2003). Stellar radiative forces will 
mostly act on the surrounding gas and dust mixture in the first 
transition region from optically thin (e.g. where the dust is evap- 
orated) to optically thick (e.g. a massive accretion disk). This 
is exactly the region, where the FLD approximation is incor- 
rect (e.g. Boley et al. 2007). To avoid this disadvantage of the 
FLD approximation, we first calculate the stellar radiative flux 
through the environment including its absorption in a corre- 
sponding first order ray-tracing routine. First order means that 
the spatial distribution of the radiative flux from stellar irradi- 
ation is calculated according to its frequency dependent optical 
depth, but re-emission of photons is shifted to a gray FLD solver. 
Sources for the thermal dust emission are the absorbed energy 
from the prior stellar irradiation step and potentially additional 
heating from the hydrodynamics, e.g. due to compression of the 
gas, accretion luminosity of sink cells, or viscous heating. In 
other words, this means that instead of solving the whole radia- 
tion transfer problem either in the FLD approximation or with a 
ray- tracing technique, we just extract the first (most important) 
absorption event of the stellar irradiation from the FLD solver 
and calculate the appropriate flux in an accurate ray-tracing man- 
ner. This splitting method allows us to consider the frequency 
dependence of the stellar spectrum in a very efficient manner. 

In the first Subsect. 2. 1 , we recapitulate the FLD equation for 
thermal dust emission. In the following two Subsects. 2.2 and 
2.3, we explain how this FLD solver can simply be combined 
with a first order ray-tracing routine (either gray or frequency 
dependent respectively) to include irradiation feedback from a 
single central object. Due to the fact that these rays are aligned 
with the radial coordinate axis in spherical coordinates, this kind 
of coordinate system is highly favored, but not required, to solve 
the ray-tracing step. In the last Subsect. 2.4, we comment on the 
so-called generalized minimal residual method (GMRES), our 
default implicit solver algorithm for the FLD equation. 

2.1. Flux Limited Diffusion 

The thermal dust emission is solved in the FLD approximation 
based on a diffusion equation for the thermal radiation energy 
E R . Within a given spatial density p(x) and temperature T(x) 
distribution, we start from the time evolution of the internal en- 
ergy density Ei nt and thermal radiation energy density E R : 

d,E mt + V ■ (£ int u) = -PVu+A (1) 
d,E R + V -(E R u) = -V-F-A (2) 

with the corresponding thermal pressure P, dynamical velocity 
h, and flux of radiation energy density F. The radiative heating 
and cooling of the gas is covered in A = p c k r (aT 4 - £r), 
where c is the speed of light, k r the Rosseland mean opacity, and 
a the radiation constant. Gas and dust temperatures are assumed 
to be the same (r gas = 7d us t = TJ. 
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We add up the Eqs. (1) and (2) and solve the transport term 
V {E R h) separately, if necessary, in the dynamical problem. If 
we consider the transport of the internal energy V ■ (Ei nt u) al- 
ready during the corresponding hydrodynamics step (operator 
splitting), the remaining terms yield 

d,(E mt + E R ) = -V-F + Q + , (3) 

where the source term Q + = —P V • u + . . . depends on the 
physics included and can additionally contain additional source 
terms such as accretion luminosity from sink cells or viscous 
heating. 

In the following, we use the assumption that the gas 
and radiation temperature are in equilibrium (also called one- 
temperature radiation transport), which is e.g. a widely used ap- 
proach in radiation hydrodynamics simulations of circumstellar 
disks. In a purely FLD radiation transport method, this assump- 
tion is justified for optically thick regions (e.g. deeply inside an 
accretion disk). The usage of our split radiation scheme guaran- 
tees also the correct gas temperature in regions of dominating 
stellar irradiation (like an optically thin envelope or disk atmo- 
sphere). In residual regions, shielded from the irradiation behind 
the optically thick circumstellar disk, the gas and radiation tem- 
perature will be the same in systems where the flow time is long 
compared to the gas-radiation equilibrium time only. The flow 
time is fg ow w l/u, where I and u are the characteristic length and 
velocity scales of the problem. The gas-energy equilibrium time 
is of the order of f eq » (a: p c) _1 , since the energy exchange rate 
is E R « k p c E R . Thus, the requirement that ffl ow » f eq reduces 
to the requirement that t » u/c, where t ~ k p I is the optical 
depth of the system. This is the regime that Mihalas & Mihalas 
(1984) identifies as static diffusion. Thus, the scheme is viable in 
the static diffusion limit. This is consistent with our test results 
for the radiative shock problem presented in Sect. 5. In this spe- 
cial case of a strong shock in a shielded, optically thin region, the 
one-temperature approach yields the correct temperature of the 
shocked gas, but results in a steeper temperature gradient than a 
two-temperature radiation transport method in the upstream di- 
rection of the radiation flux. Consideration of such small effects 
at outer regions of the domain goes beyond the scope of this re- 
search project, which focuses on the details of stellar radiative 
feedback on the accretion flow onto a massive star. Additionally, 
Krumholz et al. (2007b) show that massive star formation is gen- 
erally in the static diffusion limit. In the end, we will benefit from 
the speedup of the radiation transport solver by cutting the num- 
bers of unknown variables in half. Moreover, the stiffness of the 
set of equations is much less if the above local equilibrium as- 
sumption is made. 

Expressing both energies on the left hand side of Eq. (3) in 
terms of temperature allows us to derive a relation between the 
time derivatives of the energies. The radiation energy density in 
absence of irradiation (see following Subsect. 2.2 for the case 
including irradiation) is given by 

E R = a T\ (4) 

This expression plus the internal energy density E m — cy p T 
with the specific heat capacity cy yield 

d,E int = c y p d,T = Cy — = d,E R . (5) 
4 a T J 

With this relation, the Eq. (3) reduces to a standard diffusion 
equation 

d,E R = M-V ■ F + Q + ) (6) 
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with f c - (cy pi (A a r 3 ) + l) , depending only on the ratio 
of internal to radiation energy. The flux F of radiation energy 
density is determined in the FLD approximation via 



F = -D V£ R 



(7) 



with the diffusion constant D — Ac/ (k r p). The flux limiter A is 
chosen according to Levermore & Pomraning (1981). Scattering 
is neglected. In the most extreme limits the flux becomes either 
F — -c E R V£'r/|V£ , r| for highly optically thin regions (free- 
streaming limit) or F — -c V£r/ (3 k r p) for highly optically 
thick regions (diffusion limit). 

Summing up, we can describe the thermal/diffuse radiative 
processes via: 



d,£R = -/c(V-F-£ + ) 



(8) 



with 



fc = 


(c v pi (4 a T 3 ) + l) 


(9) 


F = 


-D V£ R 


(10) 


D = 


A c/ (k r p) 


(11) 


A = 


(2+R)/(6 + 3R + R 2 ) 


(12) 


R = 


\VE r \/(k rP E r ) 


(13) 


Q + = 


-P V ■ u + additional source terms 


(14) 



2.2. Irradiation 

When including stellar irradiation, some of the equations have 
to be modified. The irradiation from a central object is treated 
as an additional flux F», released along rays in the radially out- 
ward direction. E R is the energy density of the radiation emitted 
by the circumstellar material, and does not include the stellar 
irradiation. The additional radiation power -V ■ F« from this ir- 
radiation is therefore added to the diffusion Eq. (8) as a source 
term 



d,E R = -f c (V-F + V-F t -Q+) 



(15) 



and has to be added to the left hand side of Eq. (4), which is 
used to calculate either the thermal radiation density or the corre- 
sponding dust temperature T under the assumption that the dust 
is in equilibrium with the combined stellar and diffuse radiation 
field: 



af 



Kp(T t ) \F t \ 
k p (T) c 



(16) 



Here Kp(T) represents the Planck mean opacity to a given tem- 
perature T and T* is the star's effective temperature. 

After solving the diffusion Eq. (15) implicitly (with the ra- 
diation energy E R (x) as the solution vector), the corresponding 
dust temperature is calculated using Eq. ( 1 6). Since the opacity 
kp (T) depends on the temperature the right hand side of Eq. (16) 
also depends on the solution of T. This makes an iterative pro- 
cedure based on the Newton-Raphson method necessary to find 
the solution. The ratio of opacities a:p (T*) /k r (T) corresponds to 
the relation of emission from and absorption by dust. 

The stellar radiative flux as a function of distance r from the 
central object is calculated by 



(R \ 
y) exp(-r(r)) 



(17) 



The (boundary) flux at the stellar surface F* (R*) can be calcu- 
lated from the Stefan-Boltzmann law 



IF. (R.) | = cr SB Tt 



(18) 



with the Stefan-Boltzmann constant <tsb — 0.25 a c. The op- 
tical depth t is, in the case of gray irradiation (see following 
Subsect. 2.3 for specific changes due to frequency dependent ir- 
radiation), calculated along radial rays through the spherical grid 
via 



t(D 



f 

Jr. 



K P (T t ) p(r) dr. 



(19) 



In the case of the radiation benchmark test by Pascucci et al. 
(2004), the optical depth between the stellar surface and the in- 
ner disk (the inner boundary of the spherical computational do- 
main) is assumed to be negligible. 

A flow chart of the radiation module described so far is 
shown for a static problem such as the benchmark test of Sect. 3 
in Fig. 1. 



Stellar parameter: 

T„R, 



Density distribution: 

pW 



Initial Setup 



Irradiation: 

F",(T„R„K P (T,),p(x)) 



from Eq. (17) 



Radiation Energy: 

E£(r\F;) 



from Eq. (16) 



implicit FLD solver 



New Radiation Energy: 
£7' [df, El, k r (x)) from Ec l- ( 15 ) 

iterative Newton-Update 



Main Loop 



New Temperature distribution 

T" +l 




from Eq. ( 1 6) 



Fig. 1. Schematic flow chart of the radiation module. Exponents 
declare the iteration (timestep) number n. The actual timestep 
corresponds to df = f -f~ x . 



2.3. Frequency dependent irradiation 

Taking into account the frequency dependence of the stellar flux, 
we consider a fixed number of frequency bins, characterized by 
their mid-frequency v,, instead of the gray irradiation. For the 
radiation benchmark test we use the opacity tables of Draine & 
Lee (1984), including 61 frequency bins, see Fig. 2. Each dot 
in Fig. 2 represents the mid-frequency of the corresponding fre- 
quency bin. Including the frequency dependence, we have to re- 
place the Planck mean opacity k r (T„) in Eqs. (16) and (19) with 
the frequency dependent opacities k(v) and the stellar flux F» 
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Fig. 2. Regarding the frequency dependence of stellar irradia- 
tion feedback. Frequency dependent opacities k (v), Planck mean 
opacities k p (T t ), and stellar black body spectrum B v (v, T t ) as 
function of the frequency v. 



in Eqs. (15) and (16) with the corresponding sum over the fre- 
quency dependent fluxes F* (v), each calculated via Eq. (17). 

The boundary condition, given by Eq. (18), is now deter- 
mined for each frequency bin by the integral over the corre- 
sponding part of the stellar black body Planck function B v (v, T t ) 
in the frequency dependent interval: 



(v,-+v, +1 )/2 



I 



B y (v, r.) dv 



(20) 



(v;-l+Vi)/2 



2.4. The generalized minimal residual method 

The FLD Eqs. (8) or (15) are solved in our newly implemented 
approximate radiation transport module for hydrodynamics in an 
implicit fashion via a generalized minimal residual (GMRES) 
solver. The GMRES solver is a Krylov subspace (KSP) method 
for solving a system of linear equations Ax = b via an ap- 
proximate inversion of the large but sparse matrix A and is an 
advancement of the minimal residual (MinRes) method. The 
GMRES method was developed in 1986 and is described in Saad 
& Schultz (1986). This method is much better than the conju- 
gate gradient (CG) or successive over-relaxation (SOR) method, 
and at least as good as the well-known improved stabilized Bi- 
Conjugate Gradient (BiCGstab) method used, e.g., in Yorke & 
Sonnhalter (2002). 

The general idea of minimal residual solvers based on 
the KSP method is the following: The / th KSP is defined as 
Ki = span{b,Ab,A 2 b,...,A'~ l b}. In each solver iteration i the 
GMRES method increments the subspace K, used with an ad- 
ditional basis vector A' -1 b and approximates the solution of the 
system of linear equations by the vector Xj which minimizes the 
norm of the residual r - \Ax, - b\. This method converges mono- 
tonically and theoretically reaches the exact solution after per- 
forming as many iterations as the column number of the matrix 
A (which equals the number of grid cells). Of course, in practice 
the iteration is already stopped after reaching a specified relative 
or absolute tolerance of the residual, which normally takes only 
a small number of iterations. The computation of each iteration 
grows like O (; 2 ). In the current implementation, we use the so- 
called 'GMRES restarted' by default. GMRES restarted never 
performs all the iterations to reach the exact solution. After a 
priorily fixed number of internal iterations n, the solver starts a 



second time in the first subspace K\ but with the last approximate 
solution x„. Due to the growing computational effort with O (irj 
this approach generally results in a speedup of the computation. 

The radiation transport module is parallelized for distributed 
memory machines, using the message passing interface (MPI) 
language. The results of a detailed parallel performance test of 
the whole radiation transport module, including this GMRES 
solver is presented in Sect. 4. 

3. Frequency dependent test of the approximate 
radiation transport 

The approximate radiation transport introduced in the previous 
section can now be tested for realistic dust opacities in a stan- 
dard benchmark test for irradiated circumstellar disk models. 
The setup of the following comparison is adopted from Pascucci 
et al. (2004) and includes a central solar-type star, an irradiated 
circumstellar flared disk and an envelope. We have to choose 
a low-mass central star, because no benchmark for high-mass 
stars was performed so far. However, the tests should not depend 
on the actual size and luminosity of the central star. For com- 
parison, we choose a standard full frequency dependent Monte- 
Carlo based radiation transport code. The comparison is done 
for two different (low and high) optical depths taken from the 
original radiation benchmark test. To test each of the compo- 
nents (gray and frequency dependent irradiation as well as Flux 
Limited Diffusion) of the proposed approximate radiative trans- 
fer separately, we perform several test runs with and without the 
different physical processes (absorption and re-emission) with 
both the Monte-Carlo based (see table 1) and the approximate 
(see table 2) radiative transfer code. 

3.1. Setup 

3.1 .1 . Physical setup of the star, the disk, and the envelope 

The stellar parameters are solar-like: The effective temperature 
of the star is 5800 K (T t = 5800 K) and the stellar radius is fixed 
to 1 solar radius (R* = 1 R ). The disk ranges from r m j n = 1 AU 
up to r max = 1000 AU. 

Although the numerical setup of the gas density is done in spher- 
ical coordinates, the analytic setup of the gas density, as de- 
scribed in the original benchmark test, is given in cylindrical 
coordinates: 



P (r, z) = pa f\ (r) f 2 (r, z) 
with the radially and vertically dependent functions 

Mr) 



r 



and 



f 2 (r, z) = exp 



z 

W) 



2\ 



making use of the following abbreviations 

, 1.125 

h(r) = z d I- 
Vd 

r d = -^p = 500 AU, and 



Zd 



125 AU. 



(21) 
(22) 

(23) 

(24) 

(25) 
(26) 
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The lowest density is limited to a relative factor of 10 
compared to the highest density (at r m ; n in the midplane) to 
avoid divisions by zero (e.g. in the calculation of the diffusion 
coefficients). The normalization po of the density setup is chosen 
to define different optical depths T55o nm through the midplane of 
the corresponding circumstellar disk (at a visual wavelength of 
550 nm): 



T 550nm 


Po [g cm" 3 ] 


M tot of gas [Mq] 


0.1 
100 


8.321 * 10~ 21 
8.321 * 10~ 18 


1.1 * 10" 5 
1.1 * 10~ 2 



The opacity tables used are the same as in the original bench- 
mark test (Pascucci et al. 2004) taken from Draine & Lee (1984). 
They are displayed in Fig. 2. 

3.1.2. Numerical setup of the approximate radiation transfer 
code 

The runs of the approximate radiative transfer code are per- 
formed on a radially stretched, polar uniform, spherical, two- 
dimensional grid. The grid consists of 60 cells in the radial di- 
rection times 61 cells in the polar direction (plus additional cells 
for the storage of boundary conditions). The polar range covers 
the full spatial setup of 180°. 

Stretching of the radial grid dimension by an additional 
10% from one cell to the next is applied. The implicit diffu- 
sion Eq. (15) is solved via the GMRES method (see Sect. 2.4) 
after parallel/global Block-Jacobian and serial/local ILU pre- 
conditioning in the framework of the version 2.3.3 of the 
open source parallel solver library PETSc (Portable, Extensible 
Toolkit for Scientific Computation). More detailed information 
about this solver library can be found in Balay et al. (2001, 
2004). 

The gradient of the radiation energy is zero at the inner 
radial and both polar boundaries (V£r = 0), i.e. radiative flux 
over these boundaries is prohibited. The radial outer boundary 
is defined as a constant Dirichlet boundary corresponding to 
T = 14.7 K (E R = a T*). 

The timestep used for the FLD solver is 10 4 s. The main 
iteration (circle of irradiation and FLD steps) is stopped when 
the relative change of the temperature in each cell is smaller than 
0.01%, leading to 3 main iterations in the "purely absorption" 
runs and more than 600 main iterations in the "irradiation plus 
FLD" runs. 

3.1 .3. Numerical setup of the Monte-Carlo based 
comparison code RADMC 

For comparison we use the Monte-Carlo based radiation transfer 
code RADMC described in Dullemond & Turolla (2000) and 
Dullemond & Dominik (2004). The general solver method of 
the code is based on Bjorkman & Wood (2001). The Monte- 
Carlo runs are performed on a 60 x 3 1 grid assuming symmetry 
to the disk midplane. The grid is stretched in both directions. 
One million photons are used. 

Scattering can be handled by this code, but is simply 
switched off as it is neglected in the radiation transport module 
described here. Scattering would increase the temperature in the 
irradiated parts (up to an optical depth of about unity) by about 
2% in the optically thin envelope up to a maximum of 19% in 
the optically thick inner rim of the disks midplane due to higher 



extinction. The more effectively shielded outer regions of the 
disk would be about 10% cooler. For more massive and lumi- 
nous stars the effect of scattering would decrease due to stronger 
forward scattering, which is included in our ray-tracing routine 
per definition. 

3.2. Configurations of runs performed 

The following comparison of the results is divided into three 
parts (each of the proposed components of the theoretical 
Sects. 2.1 - 2.3 are tested independently): First, we study a pure 
absorption scenario without diffusion in the optically thin and 
thick case. Afterwards we include diffusion effects. Therefore, 
we perform three different runs of the Monte-Carlo based code: 
a full run for the optically thin case Tssonm = 0.1, a full run for 
the optically thick case T5so n m = 100, and an additional run for 
the optically thick case T5 5 o nm = 100 with excluded re-emission 
of the photons (to achieve a pure absorption scenario for com- 
parison with our first order ray-tracing routine), hereafter called 
"one-photon-limit". An overview of these three comparison runs 
is given in table 1 . 

The resulting temperature distributions of these Monte-Carlo 
runs are compared with the results of the approximate radiative 
transfer runs including different components of our module: We 
discuss five different configurations for the optically thin and 
thick setup including gray and frequency dependent irradiation 
plus potential diffusion. An overview of the physics applied and 
the resulting deviations of these runs from the comparison data 
is given in table 2. These results are discussed and illustrated in 
detail in the following subsections. 

3.3. Results 

3.3.1. Gray absorption in an optically thin disk 

In the most optically thin case Tssonm = 0.1, diffusion effects 
should be negligible. Therefore, we can test the validity of 
the routines described in Sect. 2.2 without running the diffu- 
sion routine. Compared to the full Monte-Carlo simulation from 
RADMC also the deviation in the most "difficult" region of the 
midplane (due to having the highest absorption) stays below 2% 
(see Fig. 3). 

In this optically thin limit, the diffusion effects are indeed 
negligible: An additionally performed full run of the approxi- 
mate radiation transfer module with frequency dependent irradi- 
ation plus Flux Limited Diffusion shows a variation in the radial 
temperature slope from the pure gray irradiation run below 1%. 

3.3.2. Gray and frequency dependent absorption in an 
optically thick disk 

In the most optically thick case T55o nm = 100, we test a pure 
absorption case to distinguish the deviations introduced by the 
FLD approximation in the full run from the deviations intro- 
duced by the irradiation component of our module. Therefore, 
we run the Monte-Carlo code only until every initial photon is 
absorbed, neglecting re-emission or scattering events. In this sce- 
nario, we are able to probe the absorption routines of Sects. 2.2 
and 2.3 in detail and determine the improvement by considering 
the frequency dependence of the stellar irradiation. The result- 
ing temperature profiles through the midplane for the case of 
gray and frequency dependent irradiation as well as the corre- 
sponding deviations from the Monte-Carlo comparison run are 
shown in Fig. 4. Indeed, we found that in order to limit the devi- 
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Run 


MCO.l-full 


MCI 00- A 


MClOO-full 


Optical depth r 5 5o nra 


0.1 


100 


100 


RADMC - Configuration 


"full" 


"one-photon-limit" 


"full" 


frequency dependence 


yes 


yes 


yes 


Comparison Sect. 


3.3.1 


3.3.2 


3.3.3 & 3.3.4 



Table 1. Overview of performed Monte-Carlo runs for comparison. 

The overview table of the comparison runs performed with the Monte-Carlo based code RADMC contains the corresponding optical 
depth of the test case and the used configuration of the Monte-Carlo code, see referred comparison Sects, for details on the "one- 
photon-limit" and "full run". All Monte-Carlo runs include frequency dependence. 



Run 


G0.1 


G100 


F100 


GDI 00 


FD100 


Optical depth Tssonm 


0.1 


100 


100 


100 


100 


Comparison run 


MCO.l-full 


MCI 00- A 


MCI 00- A 


MClOO-full 


MClOO-full 


Ray-tracing Irradiation 


gray 


gray 


freq. dep. 


gray 


freq. dep. 


Flux Limited Diffusion 


no 


no 


no 


yes 


yes 


Theory Sect. 


2.2 


2.2 


2.2 & 2.3 


2.1 & 2.2 


2.1 - 2.3 


Comparison Sect. 


3.3.1 


3.3.2 


3.3.2 


3.3.3 


3.3.3 & 3.3.4 


Deviation (AT) /T [%] 


< +2.0 


+ 10.9 -> -57.2 


+5.0 -> -0.6 


+5.0 -> -38.4 


+9.3 -> -11.1 



Table 2. Overview of runs using the proposed approximate radiation transport module. 

The overview table of the runs discussed contains the corresponding optical depth of the test case and the Monte-Carlo run, which is 
used for comparison (see also table 1 and referred theory and comparison Sects.). Furthermore, the applied radiative modules (gray 
or frequency dependent absorption as well as possible diffusion) of our proposed approximate radiation transport method and the 
corresponding Sects., in which this modules and the final results are discussed, are given. The deviations in the resulting temperature 
profiles of the approximate radiation transport from the corresponding Monte-Carlo comparison run are shown in the lower row. 



ation in the absorption part of the radiation module to less than 
5% it is essential to account for the frequency dependence of the 
stellar irradiation. That is due to the fact that the infrared part of 
the stellar spectrum has a lower optical depth than the UV part. 
Neglecting the frequency dependence results roughly in a 57% 
cooler disk at large radii and a 10% hotter inner rim compared 
to the Monte-Carlo based radiative transfer code. 

3.3.3. Gray and frequency dependent absorption plus 
diffusion in an optically thick disk (the complete 
problem) 

In this and the following section, we show the comparison results 
between the complete approximate radiation transport method 
with the combination of (gray or frequency dependent) irradia- 
tion and Flux Limited Diffusion and the corresponding Monte- 
Carlo simulation in the most optically thick case T 5 5o nm = 100. 

Also in this case, we found that the frequency dependent ir- 
radiation is necessary to achieve the accuracy needed for realis- 
tic radiative feedback in hydrodynamics studies (e.g. in massive 
star formation or in irradiated accretion disks). Gray irradiation 
combined with the FLD approximation leads to deviations up to 
38.4% in the resulting temperature slope (see Fig. 5). Including 
frequency dependent irradiation and FLD the radial temperature 
profiles agree within 0.7% at the inner boundary up to a maxi- 
mum deviation of 11.1% at roughly r m 200 AU. The compari- 
son between the complete radiation scheme and the correspond- 
ing Monte-Carlo run is shown in radial and polar temperature 
profiles in Figs. 5 and 6 respectively. The turnover point (at a 
polar angle from the midplane of 6 « 19°) from the optically 
thin envelope to the optically thick disk region is reproduced 
very well. 

Further remarks and discussion of these results are given in 
the last Subsect. 3.4 at the end of the comparison section. 



3.3.4. Radiative force 

Stellar radiative feedback on the dynamics of the environment 
plays a crucial role in the formation of massive stars. The heat- 
ing will probably prevent further fragmentation of the cloud 
by enhancing the Jeans mass (e.g. Krumholz et al. 2007a). 
Furthermore, the dusty environment feels the radiative force 
when absorbing the radiation due to momentum conservation, 
which potentially stops the accretion process for highly lumi- 
nous massive stars. Therefore, it is necessary to compute the 
correct radiative flux (and its derivative, the radiative force) in 
addition to the temperature distribution in such simulations. 

The radiative force density is, accordingly to Mihalas & 
Mihalas (1984), given by 

p 

fmd=PK • (27) 

c 

In our split radiation transport the total radiative flux f\ot is given 
by the sum of the irradiated stellar flux F m from the ray-tracing 
routine and the diffuse component F from the FLD solver. The 
flux calculated in the FLD approximation is F = —D V£r (see 
Sect. 2.1). 

In discretized space the opacity and density are constant over 
a single grid cell and the irradiated stellar radiative flux F, is 
calculated at the cell interfaces. To calculate the mean radiative 
force inside this cell it is necessary to integrate the radiative force 
over the cell volume (e.g. a simple ansatz of averaging only the 
stored fluxes at the interfaces towards the cell center would lead 
to unphysically high radiative forces for r » 1, i.e. for the case 
that most of the flux is absorbed on a length scale much smaller 
than the grid size). Integrating the above formula Eq. (27) over 
space (for simplicity here: an one-dimensional cartesian grid 
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Fig. 3. Radial cut through the temperature profile in the midplane 

in the most optically thin case 15500m =0.1. 

Upper panel: Radial temperature slope of the gray irradiation 

routine (solid line) and the Monte-Carlo based comparison code 

(dots). 

Lower panel: Differences between the two codes in percent. 



with an uniform grid spacing of Ax) leads to: 



i r Ax 

= —r-pK F, 
c Ax Jo 



dx 



(28) 
(29) 



The flux at position x inside the grid cell is given by the flux F l „ 
entering the cell at the interface i and the absorption of this flux 
along the length x: 

F, (x) = F[e- T 

with the optical depth r = k p x. The remaining integral yields: 



i . r AA 

/* = p k F[ I exp (— k p x) dx 

c Ax Jo 



c Ax 



Fl (exp(-^pAx)-l). 



(30) 
(31) 



Finally, the mean radiative force is therefore given by the dif- 
ference of the left and the right flux into and/or out of the cell 
respectively 



f =-l F * +1 ~ ^ 
c Ax 



(32) 
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Fig. 4. Radial cut through the temperature profile in the midplane 
in the most optically thick case T55Q nm = 100 without diffusion. 
Upper panel: Radial temperature slope of gray irradiation 
(dashed line), frequency dependent irradiation (solid line), and 
the Monte-Carlo routine in the "one-photon-limit" (dots). 
Lower panel: Deviations of the gray (dashed line) and frequency 
dependent method (solid line) from the Monte-Carlo code in per- 
cent. 



which equals in continuous space the derivative of the radiative 
flux 



/* = -- d x F t . 

c 



(33) 



Combining Eqs. (27) and (33) reflects the fact that without emis- 
sion the radiative flux is given by the differential equation 



d x F t =-pKF„. 



(34) 



In other words, the relationship in Eq. (34) shows that the origi- 
nal Eq. (27) and the derived Eq. (33) are indeed identical (in con- 
tinuous physical space), whereas the left hand side of Eq. (34) 
gives the expression easily available in discretized space. In our 
code we use a grid in spherical coordinates, thus the derivative 
of the stellar irradiation is given by 



(35) 



So the purely geometrical dilution of the flux in the radially out- 
ward direction (F.(r) cc r 2 ) does of course not contribute to the 
radiative force. 
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Fig. 5. Radial cut through the temperature profile in the midplane 
in the most optically thick case T55o n m = 100 including irradia- 
tion and Flux Limited Diffusion. 

Upper panel: Radial temperature slope of gray irradiation plus 
Flux Limited Diffusion (dashed line), frequency dependent irra- 
diation plus Flux Limited Diffusion (solid line), and the corre- 
sponding Monte-Carlo routine (dots). 

Lower panel: Deviations of the gray (dashed line) and frequency 
dependent run (solid line) from the Monte-Carlo code in percent. 



We compute the radiative force in the setup of the radiation 
benchmark test with gray and frequency dependent irradiation 
as well as the thermal radiative force for the most optically thick 
case and compare our results with the corresponding Monte- 
Carlo based run. The result is visualized in Fig. 7. The peak 
position is reproduced very well. The absolute value of the peak 
is underestimated by 3-4% only. Behind the absorption peak the 
radiative force smoothly drops down and the relative deviations 
grow radially outwards, but stay below 10%. The fraction of 
the radiative force, resulting from the thermal flux, is relatively 
small and is most important at and directly after the absorption 
peak (where most of the thermal radiation is emitted). This frac- 
tion will presumably be higher in denser environments. At larger 
radii (r > 300 AU) where the disk becomes highly optically thin 
for its own thermal radiation, the radiative force resulting from 
the thermal flux is negligible (cp. Figs. 8 and 7). 

The gray approximation and the corresponding frequency 
dependent run show only small deviations (< 5%) at the outer 
part of the disk for radii roughly larger than 200 AU (after the 
absorption of the stellar irradiation). With the setup of Pascucci 
et al. (2004) the gray approximation leads to a higher radiative 
force than the frequency dependent ones, but in general the dif- 
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Fig. 6. Polar cut through the temperature profile at r — 2 AU of 
the frequency dependent irradiation plus Flux Limited Diffusion 
run for the most optically thick case T 5 5o nm = 100. The profile 
reproduces the turnover point at a polar angle of 6 « 19° above 
the midplane from the optically thin envelope to the optically 
thick disk region. 

Solid line: Frequency dependent irradiation plus Flux Limited 
Diffusion. 

Dots: Data from the corresponding Monte-Carlo comparison 
run. 

The vertical axis covers only a small temperature range from 200 
to 260 K to better visualize the small deviations. 



ference of both methods depends strongly on the stellar luminos- 
ity. The black body spectrum of more luminous stars will shift 
to higher frequencies (cp. Fig. 2). 

Further remarks, explanations, and detailed discussion of the 
resulting radiative force and temperature profiles are given in the 
following Subsect. 3.4. 

3.4. Remarks and Analysis 

3.4.1. The Monte-Carlo comparison code 

For comparison and interpretation of the results, we should men- 
tion that in the original radiation benchmark test (Pascucci et al. 
2004) the different Monte-Carlo codes themselves differ in the 
radial temperature profile through the midplane of the circum- 
stellar disk in the most optically thick case Tsgonm = 100 by 5% 
in most of the region between 1.2 to 200 AU and up to 15% to- 
wards the outer border of the computational domain in the ra- 
dial direction. This means that the deviations in the optically 
thin case (Fig. 3) as well as the frequency dependent ray-tracing 
part of the optically thick case (Fig. 4) stay beneath the discrep- 
ancy of the different Monte-Carlo solutions. As expected, the 
direct stellar irradiation is determined highly accurately, when 
considering the frequency dependence. Therefore, the errors in- 
troduced by using the FLD approximation can be limited in the 
test throughout the irradiated regions. 

The influence of the so-called photon noise in the Monte- 
Carlo method is illustrated in the highly optically thin regions 
(\6\ > 30° from the midplane) in Fig. 6, where the temperature 
should actually be independent for large polar angles (as dis- 
played by the solid line). 
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Fig. 7. Radial cut through the radiative force profile at a polar 
angle of 8 m 27° from the midplane in the most optically thick 
case T 5 5o nrn = 100. The profile visualizes the onset of the radia- 
tive force at the transition from the optically thin envelope to the 
optically thick disk region. 

Upper panel: Radial radiative force profile of frequency depen- 
dent irradiation plus Flux Limited Diffusion (solid line with 
squares) as well as the result from the Monte-Carlo routine (solid 
line). 

Lower panel: Deviations of purely stellar (dotted line with tri- 
angles) and stellar plus thermal radiative force (solid line with 
squares) from the Monte-Carlo run in percent. 



3.4.2. The Flux Limited Diffusion approximation 

A special feature of the setup of the original benchmark test 
(Pascucci et al. 2004) is the fact that even in the most optically 
thick case Tgsonm = 100, which is defined for a wavelength of 
550 nm, the pre-described disk is locally optically thin for the ra- 
diation from thermal dust emission. Integrating the correspond- 
ing local optical depth tr (r) = kr(T) p(r) Ar from the outer 
edge of the disk through the midplane towards the center yields 
a final optical depth of tr « 0.5, see Fig. 8. 

The plot shows clearly the low optical depth for the thermal 
component of the radiation field especially in the outer part of 
the disk, which results in an underestimation of the temperature 
in the transition region at roughly r ~ 200 AU due to the over- 
estimation of the radiative flux in the outward direction by the 
FLD approximation. 

The FLD approximation is known to be valid in the most op- 
tically thin (free-streaming limit) as well as in the most optically 
thick (diffusion limit) regions only. The apparent yet surprisingly 




1000 



Fig. 8. Radial profile of the optical depth tr through the mid- 
plane in the most optically thick case Tssonm = 100 integrated 
from the outer edge of the disk inwards. 



good agreement between the Monte-Carlo based runs and the 
herein proposed radiation transport module also in the interme- 
diate region of the flared disk atmosphere (see Fig. 6) is due to 
the newly implemented direct irradiation routine which yields 
the correct flux and depth of penetration for the different fre- 
quency bins of the stellar irradiation spectrum (see Sect. 3.3.2). 
The slight underestimation of the temperature at r » 200 AU 
(see Fig. 5) in the disk midplane in the most optically thick case 
is most likely a result of an intermediate region (transition from 
optically thick to optically thin) in the outward radial direction, 
which is shielded from the direct irradiation and is not in good 
agreement with the FLD approximation. 

3.4.3. The frequency dependence 

Approximating the frequency dependence of the stellar spectrum 
by gray (frequency averaged) Planck mean opacities results in an 
overestimation of the optical depth in the infrared part and an un- 
derestimation of the absorption in the ultraviolet part of the stel- 
lar spectrum (see Fig. 2). This leads to an under- or overestima- 
tion of the radiative force onto dust grains at the first absorption 
peak depending on the stellar luminosity due to less absorption 
of the most energetic UV photons or stronger absorption of the 
photons at the peak of the stellar black body spectrum. Secondly, 
the gray approximation results in an underestimation of the re- 
sulting temperature deeply inside the disk due to absorption of 
the infrared photons already at inner disk radii. The quantita- 
tive amount of the resulting deviations depends in general on the 
specific stellar and dust properties. Due to the steep decay of the 
stellar black body spectrum at high frequencies, the difference of 
the gray and frequency dependent radiative force turns out to be 
very small in this specific setup. On the other hand, consideration 
of the frequency dependence is essential to limit the deviations 
in the resulting temperature profile to less than 11.1% (compared 
to 38.4% for gray irradiation plus FLD, see Fig. 5). These issues 
are well illustrated in Fig. 2. The figure shows the frequency 
dependent opacities of Draine & Lee (1984), the approximated 
frequency averaged value of the Planck mean opacity regard- 
ing the stellar effective surface temperature as well as the black 
body spectrum of the central star (to visualize the amount of ra- 
diative flux which is emitted per frequency bin). Each dot in the 
figure marks the mid-frequency of the correspondingly chosen 
frequency bin. The effects on the resulting temperature profile 
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and radiative force cannot be generalized easily. They depend 
on the underlying dust model and strongly on the properties of 
the central star, which yield a shift of the peak position of the 
black body spectrum in Fig. 2 according to Wien's displacement 
law Vpeak K T. In the specific setup of Pascucci et al. (2004) the 
Planck mean opacity at the black body peak position is higher 
than the frequency dependent ones, leading to a slightly higher 
radiative force. The strong overestimation of the opacity in the 
infrared regime leads to the huge discrepancy of the gray ap- 
proximation in the radial temperature profile through the disk. 



4. Parallel performance of the approximate 
radiation transport module 

The parallelization of the radiation transport scheme and the 
GMRES solver (see Sect. 2.4 for details of how the solver works) 
are taken care of by the PETSc library (Portable, Extensible 
Toolkit for Scientific computation), see also Balay et al. (1997). 
To test the parallel speedup of the implemented radiation trans- 
port module we perform two tests with an extended version of 
the circumstellar disk setup introduced in Sect. 3.1. 

We adopt the most optically thick setup for Tssonm = 100 
and expand it to three dimensions assuming axial symmetry. All 
runs include frequency dependent irradiation and Flux Limited 
Diffusion. The tests run for 10 main iterations, which consume 
the main computational effort for the approximate solver (later 
on, near equilibrium, the internal iterations needed decrease 
strongly). The number of internal iterations of the approximate 
implicit solver is fixed to 100 to guarantee the same amount of 
computation in all runs during this benchmark test. Due to the 
parallel Block-Jacobian pre-conditioner, the number of needed 
internal iterations (for a specified accuracy) normally increases 
with increasing number of processors. The precise value for the 
increase is strongly problem dependent (B. F. Smith, member of 
the development team of the PETSc library, private communica- 
tion). 

The parallel domain (the linear system of equations) is split 
in the azimuthal and polar direction only, which insures good 
speedup and efficiency. Decomposing the domain in the radial 
direction would decrease the parallel performance due to the 
fact that in the ray-tracing routine it is necessary to compute and 
therefore communicate the flux from the central sink cell to the 
outer boundary from the inside outward. Since the knowledge of 
the flux at the inner cell interface is needed to compute the flux 
at the outer interface, this method is hardly parallelizable as a 
domain decomposition. 

The measured times ti to t„ (n is the number of processors) 
represent the wall clock time per main iteration per grid cell 
without the non-recurring initialization and finalization of the 
code. We perform runs with 2 up to 64 processors due to the 
fact that single job submission is not available on the cluster we 
use. Cases, in which the local cache size would exceed the paral- 
lel decomposed problem size, are not taken into account, i.e. no 
misleading super linear speedup for high number of used pro- 
cessors is shown here. The speedup S is calculated as the ratio 
of the 'serial' run time compared to the wall clock time used 
by the parallel run: S = tilt n . The efficiency E is determined 
via E = tzl(t n n) = S jn. Each run for a specific grid and a 
specific number of processors is performed three times and av- 
eraged afterwards, but the differences of the resulting run times 
are negligible. 

All tests are performed on a 64-Bit Opteron cluster consist- 
ing of 80 nodes with two CPUs each. 
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4.1. The constant grid test 

The runs during this test are performed on a grid consist- 
ing of 64 x 64 x 256 grid cells. Each processor covers there- 
fore a (64 x 64 x 256) /« subdomain, depending on the number 
of processors used. This test shows the speedup one can gain 
when running a fixed problem on more and more processors. 
Therefore, the parallel efficiency declines stronger than in the 
following growing grid test due to the fact that with the usage 
of an increasing number of processors one lower the amount of 
computation and increase the amount of communication per sin- 
gle CPU. This means that the granularity (ratio of computation 
to communication) of the parallel problem drops strongly with 
an increasing number of processors. The resulting speedup fac- 
tors and efficiencies of this test are shown in Fig. 9. 
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Fig. 9. Speedup factors S = tilt n (upper panel) and correspond- 
ing efficiencies E = til (t n n) = S jn (lower panel) for a fixed 
(squares) and a growing problem size (stars). 



4.2. The growing grid test 

The runs during this test are performed on a grid consisting of 
32,768*n grid cells (e.g. a 64 x 64 x 8n grid), where n is again 
the number of processors used. Each processor covers therefore 
a subdomain containing 32,768 grid cells respectively during all 
runs, independent of the number of processors used. This test 
shows the efficiency achievable when using more processors to 
run a bigger problem. This is a realistic setup for our three- 
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dimensional radiative hydrodynamics studies of collapsing mas- 
sive pre-stellar cores, which are at the current limit of our clus- 
ters. The resulting speedup factors and efficiencies are shown in 
Fig. 9. 

4.3. Parallel performance results 

The constant grid test shows a clear speedup for the fixed prob- 
lem, so it is simply possible to compute a fixed problem faster by 
using more processors. The growing grid case shows a high effi- 
ciency of more than 95% during all runs. This speedup seems 
to be higher than the speedup of any actual freely available 
three-dimensional hydrodynamics code for spherical grids (well 
known in the astrophysical community, e.g. Pluto2, Flash2.5, 
and Zeus-MP1.5), which we tested during the summer 2007 in 
an accretion disk setup in hydrostatic equilibrium. In this sense, 
the resulting parallel speedup is high enough for the planned in- 
tegration of this radiation transfer module into a (magneto-) hy- 
drodynamics framework. 

Finally, we can recommend the usage of such an implicit 
modern KSP solver method leading to fast convergence of the ra- 
diative diffusion problem (at least for the setup discussed) while 
simultaneously offering high parallel efficiency. Admittedly, at- 
tention should be paid to the general fact that such an approxi- 
mate solver method strongly depends on the physical problem at 
hand as well as on the specified accuracy or abort criterion. 



5. Radiative hydrodynamics shock tests 

To test the approximate radiation transport scheme also in dy- 
namical interaction with a streaming fluid, we perform two stan- 
dard radiative shock tests. We adopt the setup of the supercriti- 
cal and subcritical radiative shock tests for the VISPHOT code 
in Ensman (1994). These radiative shock tests were already re- 
peated in tests of the TITAN code (Sincell et al. 1999a,b), ZEUS- 
2D (Turner & Stone 2001), as well as ZEUS-MP (Hayes et al. 
2006). Analytic approximations for this kind of problem are 
given by Zel'Dovich & Raizer (1967). 

To follow the motion of the gas, we solve (additionally to the 
radiation transport equations of Sect. 2) the equations of com- 
pressible hydrodynamics in conservative form 



d,p + V ■ (p u) 
d t {p u) + V (p u u) 
d,E + V ■ ((£ + P) u) 





-V P + pa 
pu ■ a 



(36) 
(37) 
(38) 



with the acceleration source term from the radiation transport 



V ■ F t D V£ R 

a = e r - k r , 

p c c 



(39) 



derived in Sect. 3.3.4. The evolution of the gas density p, the 
gas velocity u, the thermal pressure P, and the total energy E is 
computed using the open source magneto-hydrodynamics code 
Pluto3 (Mignone et al. 2007). Pluto is a high-order Godunov 
solver code, i.e. it uses a shock capturing Riemann solver within 
a conservative finite volume scheme. Our default configuration 
consists further of a Harten-Lax-VanLeer (hll) Riemann solver 
and a so-called "minmod" flux limiter using piecewise linear 
interpolation (plm) and a Runge-Kutta 2 (RK2) time integra- 
tion (predictor-corrector-method, compare e.g. van Leer 1979). 
Therefore the total difference scheme is 2 nd order accurate in 
time and space. 



The test setup describes a piston moving with supersonic 
velocity through an initially uniform, cold gas. The one- 
dimensional domain of the grid associated covers a distance of 
length lo = 7 * 10 10 cm. The iso-density in the domain is fixed to 
po = 7.78*10 -10 g crrT 3 . For testing purposes, the gas is set to be 
completely ionized, thus the mean molecular weight is p = 0.5. 
An ideal equation of state is used with y = 5/3. The opacity is 
fixed to a constant value of k — 0.4 cm 2 g _1 . The initial tem- 
perature drops linearly from 85 K at the starting position of the 
piston to 10 K at the outer boundary. The velocity uq > c s of 
the piston is used to determine the strength of the shock. While 
the piston moves through the domain the radiative energy from 
the shocked gas will stream upwards leading to a pre-heating, a 
pre-acceleration, as well as a pre-compression of the gas directly 
in front of the shock. If the temperature in this pre-heated region 
stays below the temperature of the shocked gas, it is called a 
subcritical radiative shock. If the temperature in the pre-heated 
region equals the temperature of the shocked gas, it is called a 
supercritical radiative shock. The smallest piston velocity lead- 
ing to a supercritical shock defines the critical velocity u c . 

Ensman (1994) used a Lagrangian grid moving with the pis- 
ton velocity. This setup is translated into an Eulerian grid by set- 
ting the initial velocity in the whole domain as well as the perma- 
nent velocity at the outer boundary to the negative of the piston 
velocity uq, compressing the gas at the inner reflective bound- 
ary, which represents the moving piston. In the visualization of 
our results (see Figs. 10 and 1 1) the spatial axes are retranslated 
into the non-moving frame used in the visualization by Ensman 
(1994) for the sake of comparison. The spherical coordinate sys- 
tem used at large radii (to achieve a planar geometry) by Ensman 
(1994) is translated into cartesian coordinates. We use 512 uni- 
form grid cells to cover the spatial extent of the grid. These grid 
adjustments were also used in the test of the ZEUS-MP code 
(Hayes et al. 2006). 

In this radiative shock test setup we are able to check the 
dynamical behaviour of the radiation module and the hydrody- 
namics. On the other hand, it implies a much easier treatment 
of radiation transport (due to the fact that the optical depth in 
front of the shock is practically constant) than the prior static but 
frequency dependent benchmark test of Pascucci et al. (2004). 
Therefore, only the FLD routine (as described in Sect. 2.1) is 
needed to perform the problem runs. 

We study this shock scenario in purely adiabatic as well as 
radiative hydrodynamics simulations. 

5.1. Radiative supercritical shock 

Higher velocity of the piston results in a stronger pre-heating of 
the gas directly in front of the shocked gas. When a piston veloc- 
ity higher than the critical velocity u c is used, the peak tempera- 
ture of the pre-heated region (which is equal to the temperature 
of the shocked gas) will not increase any more, but the enlarge- 
ment of the pre-heated region will extend further out. The gas 
temperature of the pre-heated region will presumably be in equi- 
librium with the radiation temperature for the most part of the 
domain. At one point of the sloping tail of the temperature dis- 
tribution in the upstream direction both temperatures will differ. 
In an equilibrium diffusion (also called one -temperature diffu- 
sion) code like ours the temperature will decline at this point 
sharply (see also Fig. 15 in Ensman (1994), which compares an 
adiabatic, an equilibrium, a non-equilibrium and a full radiation 
transport method for this setup). The reason for this sharp de- 
cline is (according to Ensman (1994)) that the radiation can only 
penetrate through the cold environment, if it heats this environ- 
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Fig. 10. Radiative supercritical shock. Resulting density, pressure, velocity, and temperature distributions at four different snapshots 
in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken (from left to right) at 
860 s, 4,000 s, 7,500 s, and 13,000 s after launching. Mostly horizontal lines at the lower border of the graphics refer to the initial 
setup. The snapshots at 4000 s are additionally marked by open circles for every 10 th grid cell to illustrate the resolution used. The 
spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The 
filled dots in the four profiles at 1 3,000 s represent characteristic sample points in the two-temperature approach read off Figs. 10-12 
in the original study by Ensman (1994). 



ment up to the equilibrium temperature. Once the radiation en- 
ergy is absorbed the radiative flux is zero. In a non-equilibrium 
radiation transport method the radiation energy can penetrate 
further into the environment simply by heating it to less than 
the equilibrium temperature. 

We choose the piston velocity to be the one used origi- 
nally by Ensman (1994), uq = 20 km s . In some of the sub- 
sequent tests mentioned above a slightly lower piston velocity 
(uo = 16 km s -1 ) for the supercritical radiative shock was used 
by the authors. 

Fig. 10 displays the resulting density, pressure, velocity, 
and temperature distributions for four different times (same as 
Figs. 10 - 12 in Ensman (1994)) as well as the initial setup. 
The analytic limit of a maximum jump in density by four in 
adiabatic shocks is reproduced. The effects of pre-heating, pre- 
acceleration, and pre-compression are clearly visible. The posi- 
tion of the peak temperature and the sharp declines in the den- 
sity, the pressure, and the velocity at the location of the shock 
front fit the one by Ensman (1994). That means, the propagation 
speed of the shock is reproduced correctly. Although the non- 
equilibrium gas temperature spike, the maximum temperature, 



cannot be reproduced with an equilibrium radiation transport 
method (the width of this gas temperature spike would always 
be less than a mean free path across), the value of the peak tem- 
perature is in good agreement with the radiation temperature of 
the corresponding non-equilibrium run by Ensman (1994) aside 
from minor geometrical effects of the modified coordinate sys- 
tem as already discussed in Hayes et al. (2006). Furthermore, 
the equilibrium temperature distribution resembles the tempera- 
ture distribution of the non-equilibrium run by Ensman (1994) 
in regions, where the gas temperature equals the radiation tem- 
perature (which is the most part of the domain in such a su- 
percritical shock). As expected, the temperature profile calcu- 
lated by our equilibrium diffusion method sharply declines at 
the boundary to the non-equilibrium region predicted by Ensman 
(1994). The outermost data points by Ensman (1994) in the pres- 
sure, the velocity, and the temperature distributions in Fig. 10 
display the long-range tail of the pre-heated region, if using a 
non-equilibrium radiation transport method. Summing up, the 
deviations of the equilibrium temperature method from a non- 
equilibrium approach stay below 10% in the pre- and post-shock 
regions. In this dynamic diffusion shock test, the outermost tail 
of the pre-heated region, where the gas and the radiation temper- 
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Fig. 11. Radiative subcritical shock. Resulting density, pressure, velocity, and temperature distributions at five different snapshots 
in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken (from left to right) 
at 350 s, 5,400 s, 17,000 s, 28,000 s, and 38,000 s after launching. The snapshot at 350 s mostly refers to the initial setup. The 
snapshots at 5,400 s are additionally marked by open circles for every 10 th grid cell to illustrate the resolution used. The spatial axes 
are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The filled dots 
in the temperature profile at 38,000 s represent characteristic sample points in the two-temperature approach read off Fig. 8 in the 
original study by Ensman (1994). 



ature is not in equilibrium, can naturally not be observed in an 
equilibrium approach. In the context of massive star formation, 

Krumholz et al. (2007b) already show that massive star forma- 
tion is generally in the static diffusion limit. 



5.2. Radiative subcritical shock 

In a radiative subcritical shock (c s < uq < u c ) the gas in front of 
the shock region is heated to a temperature lower than the tem- 
perature of the shocked gas region only. In a non-equilibrium ra- 
diation transport method the gas temperature sharply drops down 
(after the narrow temperature spike similar to the supercritical 
shock), whereas the radiation energy penetrates the cold un- 
shocked gas region and declines more smoothly (see also Fig. 8 
in Ensman (1994)). 

An equilibrium radiation transport is expected to first ig- 
nore the narrow gas temperature spike. Secondly, it drops down 
sharply when arriving the region, where the gas and radiation 
temperature are out of equilibrium. In other words, the equilib- 
rium radiation transport resembles the decline of the gas temper- 
ature distribution. 



We choose the piston velocity to be the one used originally 
by Ensman (1994), uq = 6 km s _I . Fig. 11 displays the result- 
ing density, pressure, velocity, and temperature distributions for 
five different times (same as Fig. 8 in Ensman (1994)). Given 
the constraints of the equilibrium radiation transport discussed 
above the resulting distributions fully satisfy the predictions. 



6. Conclusions & Outlook 

We introduced a fast and accurate radiation transport module, 
feasible for prospective usage in hydrodynamics and magneto- 
hydrodynamics simulations of circumstellar environments. The 
basic physical principle of the proposed radiation transport 
scheme is the splitting of the full radiative transfer problem into 
an accurate first order ray-tracing component and an approx- 
imate treatment of second order effects (Wolfire & Cassinelli 
1986; Murray et al. 1994; Edgar & Clarke 2003). The method 
combines thereby the advantages of the Flux Limited Diffusion 
approximation (speed) and a first order frequency dependent ray- 
tracing routine (accuracy). The thermal radiation by dust grains 
is solved in the static diffusion limit via a FLD solver. The FLD 
approximation yields the correct solution in the optically thin 
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(free-streaming) and optically thick (diffusion) limit. But in the 
transition regions, where also the most important feedback due 
to stellar radiative forces occurs, the FLD approximation is not 
valid anymore. This first impact of stellar irradiation is therefore 
accurately solved in the ray-tracing routine. 

Our implementation of the FLD method uses a MPI paral- 
lelized implicit GMRES restarted solver yielding high speedup 
factors comparable to modern hydrodynamics solvers. The ray- 
tracing routine, which calculates the distribution of stellar irradi- 
ation, benefits from the usage of spherical coordinates, in which 
the outgoing stellar rays are aligned with the radial coordinate 
axis. Expanding this first order ray-tracing to frequency depen- 
dent irradiation is therefore convenient and results in higher ac- 
curacy of the temperature slope deeply inside the disk due to 
lower optical depth of the infrared part of the stellar spectrum. A 
minor flaw of the proposed radiative transfer module is the ne- 
glecting of scattering. But in the radiation benchmark test scat- 
tering has only a marginal influence on the resulting tempera- 
ture distribution: Scattering would increase the temperature in 
the irradiated parts (up to an optical depth of about unity) due to 
higher extinction from about 2% in the optically thin envelope 
up to a maximum of 19% in the optically thick inner rim of the 
disks midplane. The more effectively shielded outer regions of 
the disk would be about 10% cooler. Secondly, for more mas- 
sive and luminous stars the effect of scattering will decrease due 
to stronger forward scattering of the most energetic UV photons, 
which is automatically included in our frequency dependent first 
order ray-tracing routine. Rowan-Robinson (1980) stated that 
neglecting the effect of scattering except for the first scattering 
of the starlight yield a good approximation. Another flaw is the 
frequency averaged instead of frequency dependent re-emission 
of once absorbed stellar photons, i.e. the expansion of the ray- 
tracing routine to second and higher order would result in a more 
and more realistic (but also computationally more expensive) ra- 
diative transfer. 

We tested the approximate radiation transport scheme for 
hydrodynamics within the setup of a well established standard 
benchmark test for (Monte-Carlo and finite differences based) 
radiative transfer codes (Pascucci et al. 2004), consisting of a 
central star, a circumstellar flared disk, and an envelope. The de- 
viations of the resulting temperature slopes stay below a maxi- 
mum of 1 1 . 1 % compared to a frequency dependent Monte -Carlo 
based radiation transport method, which requires much more 
CPU time. This deviation in temperature will be reflected in the 
associated hydrodynamics in roughly 6% deviation in the corre- 
sponding sound speed c s oc y7\ Neglecting the frequency de- 
pendence of the stellar irradiation however results in significant 
higher deviations up to 38.4% in the temperature distribution. 
Also the computed (frequency dependent) radiative force onto 
dust grains is in good agreement (below 3% at the inner rim and 
below 5% throughout the first absorption peak up to 10% rel- 
ative deviation at the outer boundary, where the radiative force 
gets negligible) with the comparison radiation transport method 
likely due to the usage of direct stellar ray-tracing feedback. The 
overall resulting deviations potentially passed to the hydrody- 
namics solver are indeed smaller than the lack of knowledge of 
the appropriate opacities. Summing up, the accuracy of the pro- 
posed radiation transport method legitimates the approximations 
needed for adequate speed in combination with hydrodynamics. 

In shock simulations, we show that the one-temperature ap- 
proach, which is valid in the static diffusion limit, is sufficiently 
accurate for the astrophysical applications we have in mind. 
Krumholz et al. (2007b) already show that massive star forma- 
tion is generally in the static diffusion limit. Nevertheless an ex- 



tension to a two-temperature approach is straight forward, if de- 
sired for other astrophysical applications in the future. 

Furthermore, we performed parallel performance tests of the 
radiative transfer module with different grid and parallelization 
setups. The parallel speedup of the solver imposes no restriction 
on large scale radiative (magneto-) hydrodynamics simulations. 

Therefore, we conclude that frequency dependent radiative 
feedback in hydrodynamics studies of the formation of a mas- 
sive star becomes feasible with the herein proposed approx- 
imate radiation transport method. Hydrodynamics simulations 
of monolithic pre-stellar core collapse models in massive star 
formation, including direct frequency dependent stellar feed- 
back, is now one of our primary goals for this radiation trans- 
port module. Secondly, three-dimensional ideal and non-ideal 
MHD simulations of the magneto-rotational instability (MRI) 
in proto-planetary accretion disks with radiative cooling are in- 
tended (Kuiper et al. 2009). Thirdly, the formation of circum- 
stellar disks and the evolution of their morphology in terms 
of potential planet formation can be studied also around low- 
mass stars. The algorithm is fully implemented in our version 
of the magneto-hydrodynamics code Pluto (including high-order 
Godunov methods for hydrodynamics, full implementation of 
the viscosity stress tensor, and additionally external forces, see 
Mignone et al. (2007)). 
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